/*
 * Ass1Prob1.c
 *
 *  Created on: 2010-4-6
 *      Author: Chunyan
 */

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#include"clcg4.h"

#define NUM_STREAMS 1  /* Number of random streams needed */
#define Z 1.96         /* alpha = 0.05 */

typedef enum {PartA, PartB, PartC,PartD,PartE} EXECUTE_MOD;

void Function_A ();
void Function_B ();
void Function_C ();
void Function_D ();
void Function_E ();
void Gen_delta2 (int n, double *pdelta_n);
void Gen_delta1 (int n, double *pdelta_n);
void Gen_delta1min2 (int n, double *pdelta_n);
void Gen_delta1div2 (int n, double *pdelta_n);
int Gen_N ();
void Gen_sample (int *sample_sn, int sample_no, int N);
int Get_N1 (int * sample_sn, int sample_no);
int Get_N2 (int * sample_sn, int sample_no);
double Cal_mu (double * pdelta, int n);
double Cal_sd (double * pdelta, int n);

void main (int argc, char * argv[])
{
  int streem;
  EXECUTE_MOD exe_mod;
  if (argc < 1) {/* missing repetition count */
     fprintf(stderr,"\nUsage: Ass1Prob1 n [-dm -t]\n");
     fflush(stderr);
     exit(0);
   }
   else {
	 if (strncmp(argv[1],"PartA",strlen(argv[1]))==0)
       exe_mod = PartA;
	 else if (strncmp (argv[1],"PartB",strlen(argv[1]))==0)
	   exe_mod = PartB;
	 else if (strncmp (argv[1],"PartC",strlen(argv[1]))==0)
	   exe_mod = PartC;
	 else if (strncmp (argv[1],"PartD",strlen(argv[1]))==0)
	   exe_mod = PartD;
	 else if (strncmp (argv[1],"PartE",strlen(argv[1]))==0)
	   exe_mod = PartE;
   }
  InitDefault ();
  for (streem = 1; streem <= NUM_STREAMS; streem++)
  {
 	InitGenerator(streem, NewSeed);
  }
 /* Correspond to different parts of the problem */
  switch (exe_mod)
  {
    case PartA :
      Function_A ();
      break;
    case PartB :
      Function_B ();
      break;
    case PartC :
      Function_C ();
      break;
    case PartD :
      Function_D ();
      break;
    case PartE :
      Function_E ();
      break;
    default :
    {
      printf ("Input Error\n");
      exit (0);
    }
  }

   exit (1);
}

/* Function for Part A  */
void Function_A ()
{
  double mu_100 = 0.0, sd_100 = 0.0; /* Expectation and variance for delta_100 */
  double mu_n = 0.0, sd_n = 0.0;     /* Expectation and variance for delta_n */
  /* Run 100 times and compute delta */
  double * pdelta_100 = malloc (sizeof (double) * 100);
  Gen_delta1 (100, pdelta_100);
//  printf ("%lf\n", pdelta_100[0]);
  int  n =0;
  double c_lower, c_upper;       /* Upper and lower level of confidence level 95% */

  /* Compute expectation and variance for delta */
    mu_100 = Cal_mu (pdelta_100, 100);
    sd_100 = Cal_sd (pdelta_100, 100);
 //   printf ("%lf %lf\n", sd_100, mu_100);


  /* Calculate n required for 5% accuracy */
    n = (int) (sd_100*sd_100*Z*Z/0.05/0.05/mu_100/mu_100) + 1;
 //   printf ("%ld\n", n);

  /* Run n times and compute delta */
   double * pdelta_n = malloc (sizeof (double) * n);
   Gen_delta1 (n, pdelta_n);
   mu_n = Cal_mu (pdelta_n, n);
   sd_n = Cal_sd (pdelta_n, n);
 //  printf ("%lf %lf\n", sd_n, mu_n);
  c_lower = mu_n - Z * sd_n/ sqrt((double)n);
  c_upper = mu_n + Z * sd_n/ sqrt((double)n);

  free (pdelta_100);
  free (pdelta_n);

  printf ("For method 1:\n");
  printf ("Required n to reach 5 percent accuracy is: %ld \n", n);
  printf ("Under n replications, 95 confidence level gives [%lf, %lf].\n", c_lower,c_upper);

}

/* Function for Part B */
void Function_B ()
{
	 double mu_100 = 0.0, sd_100 = 0.0; /* Expectation and variance for delta_100 */
	  double mu_n = 0.0, sd_n = 0.0;     /* Expectation and variance for delta_n */
	  /* Run 100 times and compute delta */
	  double * pdelta_100 = malloc (sizeof (double) * 100);
	  Gen_delta2 (100, pdelta_100);
//	  printf ("%lf\n", pdelta_100[0]);
	  int  n =0;
	  double c_lower, c_upper;       /* Upper and lower level of confidence level 95% */

	  /* Compute expectation and variance for delta */
	    mu_100 = Cal_mu (pdelta_100, 100);
	    sd_100 = Cal_sd (pdelta_100, 100);
//	    printf ("%lf %lf\n", sd_100, mu_100);


	  /* Calculate n required for 5% accuracy */
	    n = (int) (sd_100*sd_100*Z*Z/0.05/0.05/mu_100/mu_100) + 1;
//	    printf ("%d\n", n);

	  /* Run n times and compute delta */
	   double * pdelta_n = malloc (sizeof (double) * n);
	   Gen_delta2 (n, pdelta_n);
	   mu_n = Cal_mu (pdelta_n, n);
	   sd_n = Cal_sd (pdelta_n, n);
	 //  printf ("%lf %lf\n", sd_n, mu_n);
	  c_lower = mu_n - Z * sd_n/ sqrt((double)n);
	  c_upper = mu_n + Z * sd_n/ sqrt((double)n);

	  free (pdelta_100);
	  free (pdelta_n);

	  printf ("For method 2:\n");
	  printf ("Required n to reach 5 percent accuracy is: %d\n", n);
	  printf ("Under n replications, 95 confidence level gives [%lf, %lf].\n", c_lower,c_upper);

}

/* Function for Part C */
void Function_C ()
{
	// printf ("Start\n");
	 int n = 1000;  /* Iteration count */
	 double mu1_1000 = 0.0, sd1_1000 = 0.0;   /* Expectation and variance for delta1_1000 */
	 double mu2_1000 = 0.0, sd2_1000 = 0.0;   /* Expectation and variance for delta2_1000 */
	 double * pdelta1_1000 = malloc (sizeof (double) * 1000);    /* Storage for delta1_1000 */
	 double * pdelta2_1000 = malloc (sizeof (double) * 1000);    /* Storage for delta2_1000 */
	 double c1_lower, c1_upper;         /* Confidence interval for method 1 */
	 double c2_lower, c2_upper;         /* Confidence interval for method 2 */
	 double c12_lower, c12_upper;       /* Confidence interval for method mu1 minus mu2 */

	 /* Generate and calculate confidence interval of delta1 */
	 Gen_delta1 (n, pdelta1_1000);
	 Gen_delta1 (n, pdelta1_1000);
     mu1_1000 = Cal_mu (pdelta1_1000, n);
     sd1_1000 = Cal_sd (pdelta1_1000, n);
	 c1_lower = mu1_1000 - Z * sd1_1000/ sqrt((double)n);
	 c1_upper = mu1_1000 + Z * sd1_1000/ sqrt((double)n);

	 /* Use new seed */
	 InitGenerator(1, NewSeed);

	 /* Generate and calculate confidence interval of delta1 */
	 Gen_delta2 (n, pdelta2_1000);
	 Gen_delta2 (n, pdelta2_1000);
	 mu2_1000 = Cal_mu (pdelta2_1000, n);
	 sd2_1000 = Cal_sd (pdelta2_1000, n);
	 c2_lower = mu2_1000 - Z * sd2_1000/ sqrt((double)n);
	 c2_upper = mu2_1000 + Z * sd2_1000/ sqrt((double)n);

	 printf ("For method 1:");
	 printf ("Under 1000 replications, 95 confidence level gives [%lf, %lf].\n", c1_lower,c1_upper);
	 printf ("For method 2:");
	 printf ("Under 1000 replications, 95 confidence level gives [%lf, %lf]\n", c2_lower,c2_upper);

	 /* Calculation for method 1 minus method 2 */
	 c12_lower = mu1_1000 - mu2_1000 - Z * sqrt(sd1_1000*sd1_1000 + sd2_1000*sd2_1000)/ sqrt((double)n);
	 c12_upper = mu1_1000 - mu2_1000 + Z * sqrt(sd1_1000*sd1_1000 + sd2_1000*sd2_1000)/ sqrt((double)n);

	 printf ("For method 1 minus method 2:");
	 printf ("Under 1000 replications, 95 confidence level gives [%lf, %lf]\n", c12_lower,c12_upper);

	 free (pdelta1_1000);
	 free (pdelta2_1000);
}


/* Function for Part D */
void Function_D ()
{
  double mu12_1000 = 0.0, sd12_1000 = 0.0; /* Expectation and variance for delta12_1000 */
  double * pdelta12_1000 = malloc (sizeof (double) * 1000);
  double c12_lower, c12_upper;
  int n = 1000;

  /* Generate 1000 difference and calculate */
  Gen_delta1min2 (n,pdelta12_1000);
  mu12_1000 = Cal_mu (pdelta12_1000, n);
  sd12_1000 = Cal_sd (pdelta12_1000, n);
  c12_lower = mu12_1000 - Z * sd12_1000/ sqrt((double)n);
  c12_upper = mu12_1000 + Z * sd12_1000/ sqrt((double)n);

  printf ("For the difference of expected errors:\n");
  printf ("Under %d replications, 95 confidence level gives [%lf, %lf]\n", n, c12_lower,c12_upper);

  free (pdelta12_1000);

}

/* Function for Part E */
void Function_E ()
{

  double mu12_2000 = 0.0, sd12_2000 = 0.0; /* Expectation and variance for delta12_1000 */
  double * pdelta12_2000 = malloc (sizeof (double) * 2000);
  double c12_lower, c12_upper;
  int n = 2000;

  /* Generate 1000 difference and calculate */
  Gen_delta1div2 (n,pdelta12_2000);
  mu12_2000 = Cal_mu (pdelta12_2000, n);
  sd12_2000 = Cal_sd (pdelta12_2000, n);
  c12_lower = mu12_2000 - Z * sd12_2000/ sqrt((double)n);
  c12_upper = mu12_2000 + Z * sd12_2000/ sqrt((double)n);

  printf ("For the ratio of expected errors:\n");
  printf ("Under %d replications, 95 confidence level gives [%lf, %lf]\n", n, c12_lower,c12_upper);

  free (pdelta12_2000);

}



/* Generate number N uniformly distributed between 10,000 and 100,000 */
int Gen_N ()
{
  return ((int)(10000 + 90000 * GenVal (1)));
}

/* Generate sample_no serial numbers from N serial numbers */
void Gen_sample (int *sample_sn, int sample_no, int N)
{
  int i, j, a ,b;
  j = 0;
  a = sample_no;
  b = N;
  double p = (double)a/(double)b;
  for (i = 0; i < N; i++)
  {
	  if (GenVal(1) < p)
	  {
	    sample_sn[j] = i;
	//    printf ("%d\n", i);
	    j++;
	    a -= 1;
	    if (a == 0)
	    	break;
	  }
	  b -= 1;
	  p = (double)a/(double)b;

  }
}

/* Generate n Delta using method 1*/
void Gen_delta1 (int n, double *pdelta_n)
{
  int sample_no, N, N1, i;
  for (i = 0; i < n; i++)
  {
	  N = Gen_N ();
	  sample_no = (int)(0.05 * N) + 1;
//	  printf ("%ld %ld\n", N, sample_no);
	  int *sample_sn = malloc (sizeof (int) * sample_no);
	  Gen_sample (sample_sn, sample_no, N);
	  N1 = Get_N1 (sample_sn, sample_no);
	  if (N1 > N)
	    pdelta_n[i] = (double)(N1 - N)/N*100.0;
	  else
		pdelta_n[i] = (double)(N - N1)/N*100.0;
	//  printf ("%lf\n",pdelta_n[i]);
	//  	  printf ("%ld %ld\n", N, N1);
	  free (sample_sn);

 }
 // return (pdelta_n);
}

/* Generate n Delta using method 2*/
void Gen_delta2 (int n, double *pdelta_n)
{
  int sample_no, N, N2, i;
  for (i = 0; i < n; i++)
  {
	  N = Gen_N ();
	  sample_no = (int)(0.05 * N) + 1;
//	  printf ("%ld %ld\n", N, sample_no);
	  int *sample_sn = malloc (sizeof (int) * sample_no);
	  Gen_sample (sample_sn, sample_no, N);
	  N2 = Get_N2 (sample_sn, sample_no);
	  if (N2 > N)
	    pdelta_n[i] = (double)(N2 - N)/N*100.0;
	  else
		pdelta_n[i] = (double)(N - N2)/N*100.0;
	//  printf ("%lf\n",pdelta_n[i]);
	// 	  printf ("%ld %ld\n", N, N2);
	  free (sample_sn);

 }
}

/* Generate the difference between delta1 and delta2 */
void Gen_delta1min2 (int n, double *pdelta_n)
{
  int sample_no, N, N1, N2, i;
  for (i = 0; i < n; i++)
  {
	  N = Gen_N ();
	  sample_no = (int)(0.05 * N) + 1;
	  double temp1, temp2;
	  int *sample_sn = malloc (sizeof (int) * sample_no);

	  Gen_sample (sample_sn, sample_no, N);
	  N1 = Get_N1 (sample_sn, sample_no);
	  N2 = Get_N2 (sample_sn, sample_no);

	  if (N1 > N)
	  	  temp1 = (double)(N1 - N)/N*100.0;
	  else
	  	  temp1 = (double)(N - N1)/N*100.0;

	  if (N2 > N)
	      temp2 = (double)(N2 - N)/N*100.0;
	  else
		  temp2 = (double)(N - N2)/N*100.0;

	  pdelta_n[i] = temp1 - temp2;

	  free (sample_sn);

 }
}

/* Generate the ratio between delta1 and delta2 */
void Gen_delta1div2 (int n, double *pdelta_n)
{
  int sample_no, N, N1, N2, i;
  for (i = 0; i < n; i++)
  {
	  N = Gen_N ();
	  sample_no = (int)(0.05 * N) + 1;
	  double temp1, temp2;
	  int *sample_sn = malloc (sizeof (int) * sample_no);

	  Gen_sample (sample_sn, sample_no, N);
	  N1 = Get_N1 (sample_sn, sample_no);
	  N2 = Get_N2 (sample_sn, sample_no);


	  if (N1 > N)
	  	  temp1 = (double)(N1 - N)/N*100.0;
	  else
	  	  temp1 = (double)(N - N1)/N*100.0;

	  if (N2 > N)
	      temp2 = (double)(N2 - N)/N*100.0;
	  else
		  temp2 = (double)(N - N2)/N*100.0;

	  if (temp2 != 0)
	     pdelta_n[i] = temp1 / temp2;
	  else
		 i--;
//	  printf ("%lf %lf %lf\n",temp1, temp2, pdelta_n[i]);
	  free (sample_sn);

 }


}


/* Get N1 from serial number sample sample_sn */
int Get_N1 (int * sample_sn, int sample_no)
{
    int m = 0;
	if (sample_no%2 == 0)
	{
		m = sample_no/2 - 1;
	//	printf ("%d %d\n", sample_no, m);
		return (sample_sn[m]+sample_sn[m+1]-1);
	}
	else
	{
		m = (int)(sample_no/2);
	//	printf ("%d %d\n", sample_no, m);
		return (2*sample_sn[m] - 1);
	}

}

/* Get N2 from serial number sample sample_sn */
int Get_N2 (int * sample_sn, int sample_no)
{
   int m;
   m = sample_sn[sample_no-1];
   float temp = (sample_no + 1)/sample_no * m - 1;
   return ((temp-(int)temp)>0.5 ? (int)temp+1:(int)temp);
}

/* Calculate expectation of delta */
double Cal_mu (double * pdelta, int n)
{
  double mu = 0.0;
  int i;
  for (i = 0; i < n; i++)
  {
	  mu += pdelta[i];

  }

  return (mu/(double)n);
}

/* Calculate standard deviation of delta */
double Cal_sd (double * pdelta, int n)
{
	double sd = 0.0;
	double mu = 0.0;
	int i;
	mu = Cal_mu (pdelta, n);
	for (i = 0; i<n ; i++)
	{
	   sd += (pdelta[i] - mu)*(pdelta[i]-mu);
	}
  //  printf ("%lf\n", sd);
	return sqrt (sd/(double)(n-1));
}

